Force Chains, Microelasticity and Macroelasticity 
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It has been claimed that quasistatic granular materials, as well as nanoscale materials, exhibit 
departures from elasticity even at small loadings. It is demonstrated, using 2D and 3D models with 
interparticle harmonic interactions, that such departures are expected at small scales [below 0(100) 
particle diameters], at which continuum elasticity is invalid, and vanish at large scales. The models 
exhibit force chains on small scales, and force and stress distributions which agree with experimental 
findings. Effects of anisotropy, disorder and boundary conditions are discussed as well. 

PACS numbers: 45.70.Cc, 46.25.Cc, 83.80.Fg, 61.46.+W 
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There are at least two classes of systems whose appar- 
ent departure from elastic mechanical response, even for 
infinitesimal loads, has been recently discussed in the lit- 
erature: granular (lj and nanoscale materials. These 
two classes share a common property: both are typically 
composed of a relatively small number of constituents, 
macroscopic grains in the former and atoms in the latter 
case. Nanoscale grains are not considered here. The dis- 
cussion below, though relevant to both classes, is worded 
in "granular terms" and focuses on granular materials. 

Models whereby forces in static granular matter "prop- 
agate" (corresponding to hyperbolic or parabolic PDE's) 
have been suggested ||, || . This mechanism is in marked 
contrast with the non-propagating nature of the classi- 
cal equations of static elasticity (elasto-plastic models are 
successfully used by engineers in this field §,§). 

One of the goals of this Letter is to show that some 
recent experiments |7| |8j, in which the quasi-static re- 
sponse of granular matter has been measured, are con- 
sistent with an elastic description, up to differences which 
stem from the fact that continuum elasticity is a macro- 
scopic theory, valid only above a certain spatial scale. 
A crossover between microelasticity, i.e., the small scale 
response of a system whose constituents interact by har- 
monic forces, and continuum elasticity ('macroelasticity') 
at large scales is demonstrated. Note that both the re- 
sponse of small (elastic) systems and the short distance 
response of any such system are microelastic. 

Any particle in a granular medium experiences inter- 
actions with a finite number of other particles, each of 
which defines a different direction. Therefore the local 
environment of a particle is not isotropic. The existence 
of preferred directions on the particle scale implies the 
possible emergence of force chains, i.e., chains of contacts 
along which the forces are e.g., stronger than the mean 
interparticle force. This fact does not preclude elastic- 
ity, e.g., the forces may be derived from a potential (and 
can be linearized) yielding a set of discrete elliptic equa- 
tions. The notion of force propagation along force chains 
is merely an approximation pertaining to the "strong" 
forces. Continuum theory cannot be expected to de- 



scribe these microscopic interparticle forces. The macro- 
scopic field which is related to the microscopic forces is 
the stress. Even at small but finite spatial resolution, 
the stress tensor is determined by an appropriate average 
over forces, and it may bear little or no resemblance to 
the force distribution on the microscopic (particle) scale. 
Though the stress field is well defined on small scales, the 
constitutive relations correspond to continuum elasticity 
only on sufficiently large scales. In particular, disorder 
can have a large effect at small scales, yet allow for large 
scale homogeneous (and possibly isotropic) elasticity. 

The discussion below does not apply to isostatic sys- 
tems or systems at incipient failure, where particle rear- 
rangements are prominent and the range of validity of 
elasticity may be extremely limited. We consider only 
dense systems (which form the majority of granular sys- 
tems in practical applications) where geometrical con- 
straints and friction prevent major rearrangements for 
sufficiently small loads (depending on the yield stress). 

Consider a two dimensional (2D) collection of uni- 
form disks resting in a finite enclosure (with a horizontal 
'floor') under the influence of gravity. A unit downward 
vertical force is applied at the top of the system (in the 
simulations below the force acts in the middle of the up- 
per layer). Fig. |l| presents the results of a simulation in 
which the disks are assumed to interact by uniform linear 
forces ('springs' whose rest length is the particle diame- 
ter). Force chains are evident. Since non-cohesive grains 
cannot experience tensile forces, we have repeated the 
simulation with a more realistic interaction in which 'one- 
sided springs' (which can exert only compressional forces) 
have been employed. The resulting force distribution is 
depicted in Fig. |^. Though in the second simulation there 
is particle rearrangement (some contacts are severed, as 
observed in || for a pile geometry) the appearance of the 
force chains in both systems is very similar. The force 
distribution vs. the horizontal coordinate, at different 
depths, corresponding to Fig. Q, is presented in Fig. 2, 
and that corresponding to Fig. 2| is presented in Fig. 4. 
These force distributions are qualitatively similar though 
the latter (Fig. |^) is far closer to the experimental find- 
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ings H on a similar system of disks (since it better rep- 
resents the properties of the grains). Furthermore, Fig. ^ 
presents the force chains for a random network, obtained 
by adding a random number, uniformly distributed in 
the range [— |, |], to the x and y coordinates of points 
on a triangular lattice with lattice constant, d (except 
those of the bottom row), with springs connecting points 
whose distance is less than c max = 1.2d. The results are 
again qualitatively similar to those obtained in the exper- 
iment , the main difference with respect to the case of 
a regular lattice being the fact that here the force chains 
are somewhat shorter. 
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FIG. 1: Force chains in a 2D triangular lattice. A unit vertical 
force is applied to the center particle in the top layer, with no 
gravity. Line widths for all the lines shown are proportional 
to the forces. The central region of the lattice is shown (the 
entire lattice comprises 15 layers of 41 particles each). 
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FIG. 3: The norms of the interparticle forces, |f|, in the sys- 
tem depicted in Fig. |l], vs. the horizontal position, x, at sev- 
eral depths. The legend indicates the depth measured from 
the top, in layer numbers. 



0.7 
0.6 
0.5 
0.4 
0.3 



0.2 
0.1 



1 


— 




.... 2 

3 

—- 4 
6.5 


\ 




# 
■< 




* n ,* 

• 1 : 1 




1 1 f * 








* - 













x/h 



0.5 



FIG. 4: Same as Fig. H, for the case of one-sided springs. 



FIG. 2: Force chains in the same 2D triangular lattice as in 
Fig. [j], but with one-sided springs and gravity. 

Fig. ^| depicts the normal vertical stress component, 
tr zz , corresponding to the system in Fig. ^, as calculated 
using the following exact formula |l0j : 

1 f 1 

<r a p(r, t) = - ^2 fi/ja r ijp I ds(f>{r - n(t) + sr y -(i)), 








FIG. 5: Force chains in a random network of springs with 
Cmax = 1.2d (see text). 



where Ti(t) is the center of particle i, r,j(t) = rj(t)— rj(t), 
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and fj/j is the force exerted on particle i by parti- 
cle j. Greek indices denote Cartesian components. The 
coarse graining function [|l0| is <f)(r) = -ije"" 1 "!/ 1 "' , with 
w = d, the particle diameter, i.e., a fine resolution. The 
force chains are not evident any more. The model con- 
sidered here corresponds, in the continuum limit, to an 
isotropic 2D elastic medium jn| . Therefore one can com- 
pare the stress obtained from the microscopic force dis- 
tribution with that of the corresponding elastic solution. 




FIG. 6: Contour plot of ho^, in the 2D triangular lattice, 
in the region shown in Fig. pi. The contour spacing is 0.3. 
Darker shades indicate larger values of \a zz \. 

Fig. [?] compares the vertical stress at the floor of the 
system with elastic solutions for a finite slab (with rough 
or frictionless support) and a half plane. The convergence 
to the experimentally appropriate |^ (rough support) so- 
lution for a sufficient number of layers is evident. For the 
random system depicted in Fig. fi, the results are quite 
similar, expect for (expected) fluctuations. 
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FIG. 7: ha zz at the bottom of the 2D triangular lattice, com- 
pared to continuum elastic solutions. The applied force is 
unity. Note: a zz (|) = a zz (— §)• 

Consider next an anisotropic medium, obtained by tak- 
ing the spring constants for contacts in the horizontal 



direction (parallel to the slab floor), Ki, to be differ- 
ent from those for the contacts in the oblique direction, 
K 2 . As shown in Fig. ^, the obtained stress distribu- 
tion (on the floor) is either single peaked (narrower than 
the isotropic one for K 2 /Ki < 1, wider for K 2 jK\ > 1) 
or double peaked for sufficiently large K2/K1. A sim- 
ilar double peaked distribution is obtained for the case 
of 'one-sided' springs, where some horizontal contacts are 
severed, corresponding to the limit KijK\ — > 00 for these 
contacts. This is evident in Fig. |^, which clearly shows a 
macroscopic anisotropy (compare to the macroscopically 
isotropic case shown in Fig. ^). These double peaked 
distributions are similar to those obtained from hyper- 
bolic models. The results presented here indicate that 
phenomena similar to those suggested by the hyperbolic 
models can be obtained using anisotropic elasticity. The 
equations of anisotropic elasticity, which are of course 
elliptic, can approach hyperbolicity in the limit of very 
large anisotropy (as mentioned in 0). 
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FIG. 8: ha zz at the bottom of anisotropic triangular lattices 
composed of 40 layers, compared to the isotropic elastic solu- 
tion. The applied force is unity. Note: a zz (^) = a zz (— f ). 




FIG. 9: Same as Fig. g, for the case of one-sided springs. 
Very similar results to those presented above are ob- 



4 



tained for a three dimensional (3D) system: consider a 
collection of particles positioned on a simple cubic lattice. 
As shown in , this system corresponds, on large scales, 
to an isotropic elastic continuum when the spring con- 
stants for springs coupling nearest neighbors equal those 
coupling next nearest neighbors. While this model clearly 
does not describe the interaction of granular particles, we 
believe it is useful for the description of the crossover be- 
tween microelasticity and macroelasticity. Furthermore, 
as shown below, the stresses computed for this model are 
in close correspondence with the experimental results re- 
ported in |@). Consider a 3D slab of finite height, with a 
downward unit force acting on a particle in the top layer. 
As in the 2D case described above, the results for dis- 
crete lattices converge to the appropriate finite slab elas- 
tic solution with increasing depth (Fig. |l^). The effect 
of small scale anisotropy (here, the cubic symmetry) is 
seen in Fig. |ll|, which depicts contour plots of the distri- 
bution of the vertical stress on the floor. The underlying 
cubic symmetry is evident for a small depth of the system 
and is washed out for larger systems, where continuum 
elasticity holds. 
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FIG. 10: h 2 a zz along the z-axis at the bottom of a cubic 
lattice, compared to continuum elastic solutions. The results 
are scaled by h. Note: h 2 a zz (§) = h 2 a zz (-f )• 



Our results for a shallow discrete slab are similar to 
those obtained experimentally for a 'loose' granular pack- 
ing H , for which the stress distribution is narrower than 
that predicted by continuum elasticity. In contrast, for 
a dense packing, experiments show a wider distribution. 
As mentioned in Q and above, this can be explained 
by anisotropic effects (recall that the above anisotropic 
2D system can exhibit both narrower and wider distribu- 
tions). The method of preparation of the experimental 
system B suggests that a relevant model may consist of a 
number of isotropic elastic layers of variable moduli. In- 
deed, in this case the stress distribution on the floor can 
be wider or narrower than the solution for a homogeneous 
slab, depending on the distribution of moduli of the lay- 
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FIG. 11: A contour plot of h 2 a zz at the bottom of 3D slabs 
composed of 5 (solid lines) and 20 (dashed lines) discrete lay- 
ers. A unit force is applied at the top. The cubic anisotropy 
is evident for 5 layers, while the distribution appears isotropic 
for 20 layers. 



ers H In addition, the stress distribution becomes 
wider the smaller the rigidity of the support |^|, [H| . 
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